#ifndef _MCMCIMAGEFILTERS_H
#define _MCMCIMAGEFILTERS_H

#include "mcRaster.h"
#include "Vector.h"
#include "mcRasterColor.h"
#include "mcRand.h"
#include "mcmcSparseMatrix.h"
#include "mcmcHistogram.h"
#include "ItgImage.h"
#include "mcmcLowLevelStat4Affinity.h"

#define NUM_NONZERO_BINS	15

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
class MCMC_FilterBank
{
public:
	MCMC_FilterBank();
	~MCMC_FilterBank();

	virtual void	Set(Raster<float> &ras_org, const int num);
	virtual void	SetOneBank(Raster<double> &ras_filter, const int i, double *pmin=NULL, double *pmax=NULL);
	virtual void	DoFiltering(void);
	inline	int		NumFilter(void){return num_filter;};
	inline  void	NewFilterRespAtOneSite(mcmcVector<double> &vc_resp, const int k, const int j, const int i, const double value);
	inline	void	NewFilterRespAtOneSite(Histogram &histogram, const int k, const int j, const int i, const double value);
	inline	void	UpdateFilterRespAtOneSite(const int j, const int i, const double value);
	virtual void	ComputeHistograms(const int cut_off,const bool bcolor1=true);
	virtual void	FreeFilteringResult(void);
	virtual void	ComputeHistogram4Mask(mcmcMatrix<int> &mx_mask, const int k, const int cut_off, const bool bcolor1=true);
	virtual void	ExpandImage(const int row_expansion, const int col_expansion);
	virtual void	ShrinkImageBack(const int row_expansion, const int col_expansion);

protected:
	virtual void	FreeHistograms(void);

public:
	Raster<double>	*p_filtering_result;
	SparseMatrix<double> *p_filters,*p_filters_flipped;

	int				num_filter;

public:
	Raster<double>	ras_image;
	Histogram		*p_histogram_obs;
	double			*p_min_value,*p_max_value,*p_deltas;

protected:
	mcmcVector<double>		 vc_min,vc_max;
	SimpleMatrix<double> vc_predefined_minmax;
};


// filter bank responses
class MCMC_FilterResp
{
public:
	// constructions & destructions
	MCMC_FilterResp();
	~MCMC_FilterResp();
	void Set(Raster_Color<uchar> &raster1, const bool bcolor1);

	// operations
	inline int					rows(void) const;
	inline int					cols(void) const;
	inline bool					Valid(const int j, const int i) {return image_original.Valid(j,i);};
	inline	double				ImageIntensity(const int j, const int i){return ras_intensity.Data(j,i);};
	inline  double				ImageRegionLayerIntensity(const int j, const int i){return ras_region_layer.Data(j,i);};
	inline	void				GetRGB(double &R, double&G, double &B,const double in, const double p, const double t);
	__forceinline	double		ImageL(const int j, const int i){ return ras_L(j,i);};
	__forceinline	double		ImageU(const int j, const int i){ return ras_U(j, i); };
	__forceinline	double		ImageV(const int j, const int i){ return ras_V(j, i); };
	inline	double				MinL(void){return min_L;};
	inline	double				MaxL(void){return max_L;};
	inline	double				MinU(void){return min_U;};
	inline	double				MaxU(void){return max_U;};
	inline	double				MinV(void){return min_V;};
	inline	double				MaxV(void){return max_V;};
	inline  bool				GetGradient(double &g, const int j, const int i, mcmcMatrix<int> *pmx_lable=NULL);
	virtual void				SetAllFilterBanksForClustering(void);
	virtual void				SetAllFilterBanks(void);
	virtual void				ComputeHistogram(void);
	virtual void				ComputeHistogram4Mask(mcmcMatrix<int> &mx_mask, const int kcluster);
	virtual void				ComputeHistogramAsIid(RF3D &rf_iid,mcmcMatrix<int> &mx_mask, const int kcluster);
	virtual void				ComputeHistogramAsIid(RF &rf_iid,mcmcMatrix<int> &mx_mask, const int kcluster);
	virtual void				FilteringOnNewImage(Raster<uchar> &image_intensity);

			int			TotoalNumFilterResult(void);
	Raster<double>*		GetFirstResult(double &delta);
	Raster<double>*		GetNextResult(double &delta);
	virtual void		ExpandImage(const int row_expansion, const int col_expansion);
	virtual void		ShrinkImageBack(const int row_expansion, const int col_expansion);
	virtual void		DoFiltering(void);
	inline	bool		InColor(void){return bcolor;};
	inline	void		ObtainRegionResidue(irmm_Rect *pbound=NULL);
	inline  void		ObtainImageForRegion(Raster<double> &ras_curve_layer, irmm_Rect *pbound=NULL);
	virtual void		ComputeIntegralImage(void);

public:
	Raster_Color<uchar>	image_original;
	Raster<float>		ras_intensity;
	Raster<float>		ras_L, ras_U, ras_V;
	double				delta[4];
	Raster<float>		ras_filter_banks;
	Raster<float>		ras_edge_g, ras_edge_theta;
	Raster<float>		ras_DC;
	Raster<float>		ras_HF;

	MCMC_FilterBank		filter_bank_L,filter_bank_U,filter_bank_V;
	MCMC_FilterBank		filter_bank_intensity;
	int					cut_off;				// learn the parameters inside this rectangle

	// for region and curve processes
	Raster<float>		ras_region_rec;			// reconstructed image by the region layer
	Raster<float>		ras_region_residue;		// residual image by the region layer: used to adapt the
	                                            //   coefficients for each base in the curve layer
												//   = ras_intensity-ras_region_rec;
	Raster<float>		ras_region_layer;		// residual image by the curve layer: used to compute
	                                            //   the statistics of the region layer
												//   = ras_intensity-ras_region_rec;
	Raster<float>		ras_image_reconstruction;

	McItgImage			itg_image;				// integral image

protected:
	double				min_L,max_L;
	double				min_U,max_U;
	double				min_V,max_V;
	bool				bcolor;

private:
	mcmcVector<int>			vc_nbc,vc_nbr;
	int					which_result;
};

inline int	MCMC_FilterResp::rows(void) const
{
	if (bcolor)
	{
		if (filter_bank_L.ras_image.rows()>0)
			return filter_bank_L.ras_image.rows();
		else
			return image_original.rows();
	}
	else
	{
		if (filter_bank_intensity.ras_image.rows()>0)
			return filter_bank_intensity.ras_image.rows();
		else
			return image_original.rows();
	}
}

inline int	MCMC_FilterResp::cols(void) const
{
	if (bcolor)
	{
		if (filter_bank_L.ras_image.cols()>0)
			return filter_bank_L.ras_image.cols();
		else
			return image_original.cols();
	}
	else
	{
		if (filter_bank_intensity.ras_image.cols()>0)
			return filter_bank_intensity.ras_image.cols();
		else
			return image_original.cols();
	}
}

inline void MCMC_FilterResp::GetRGB(double &R, double&G, double &B, const double in, const double p, const double t)
{
	R = in*cos(p)*cos(t);
	G = in*cos(p)*sin(t);
	B = in*sin(p);
}

inline bool MCMC_FilterResp::GetGradient(double &g, const int j, const int i, mcmcMatrix<int> *pmx_label1)
{
	mcmcMatrix<int> *pmx_label;
	double		gtempL,gtempU,gtempV;
	int			l,m,n,k;

	k = 0;
	g = 0.0;
	gtempL = 0.0f;
	gtempU = 0.0f;
	gtempV = 0.0f;
	for (l=0; l<4; l++)
	{
		m = j+vc_nbr(l);
		n = i+vc_nbc(l);

		if (Valid(m,n) && (*pmx_label)(j,i)==(*pmx_label)(m,n))
		{
			gtempL += fabs((float)(image_original.L(j,i)-image_original.L(m,n)));
			gtempU += fabs((float)(image_original.U(j,i)-image_original.U(m,n)));
			gtempV += fabs((float)(image_original.V(j,i)-image_original.V(m,n)));
			k++;
		}
	}
	if (k>0)
	{
		gtempL = gtempL/((double) k);
		gtempU = gtempU/((double) k);
		gtempV = gtempV/((double) k);
		g = sqrt(gtempL*gtempL+gtempU*gtempU+gtempV*gtempV);
		if (g > 25.0)
			g = 25.0;
		return true;
	}
	else
		return false;
}

inline void MCMC_FilterResp::ObtainRegionResidue(irmm_Rect *pbound)
{
	int i,j,row_top,row_bottom,col_left,col_right;

	// get the bounding box
	if (pbound == NULL)
	{
		row_top  = 0;			row_bottom = rows()-1;
		col_left = 0;			col_right  = cols()-1;
	}
	else
	{
		row_top  = pbound->top;	row_bottom = pbound->bottom;
		col_left = pbound->left;col_right  = pbound->right;
	}
	for (j=row_top; j<=row_bottom; j++)
		for (i=col_left; i<=col_right; i++)
			ras_region_residue(j,i) = ras_intensity.Data(j,i)-ras_region_rec.Data(j,i);
}

inline void	MCMC_FilterResp::ObtainImageForRegion(Raster<double> &ras_curve_layer, irmm_Rect *pbound)
{
	int i,j,row_top,row_bottom,col_left,col_right;

	// get the bounding box
	if (pbound == NULL)
	{
		row_top  = 0;			row_bottom = rows()-1;
		col_left = 0;			col_right  = cols()-1;
	}
	else
	{
		row_top  = pbound->top;	row_bottom = pbound->bottom;
		col_left = pbound->left;col_right  = pbound->right;
	}

	double temp;
	for (j=row_top; j<=row_bottom; j++)
		for (i=col_left; i<=col_right; i++)
		{
			temp = ras_intensity.Data(j,i)-ras_curve_layer.Data(j,i);
			if (temp < INTENSITY_LOW)	temp=INTENSITY_LOW;
			if (temp > INTENSITY_HIGH)	temp = INTENSITY_HIGH;
			ras_region_layer(j,i) = temp;
		}
}


inline void MCMC_FilterBank::NewFilterRespAtOneSite(mcmcVector<double> &vc_resp, const int k, const int j, const int i, const double value)
{
	int				j1,i1,j2,i2,j3,i3;
	int				half_row,half_col,row1,row2,col1,col2;

	half_row = (p_filters_flipped[k].rows()-1)/2;
	half_col = (p_filters_flipped[k].cols()-1)/2;

	row1 = max(-half_row,0-j); row2 = min(half_row,ras_image.rows()-1-j);
	col1 = max(-half_col,0-i); col2 = min(half_col,ras_image.cols()-1-i);

	if (p_filters_flipped[k].NumSparse()<NUM_NONZERO_BINS)
	{
		// how many response changes regarding to the change at this site
		vc_resp.SetRowVector(p_filters_flipped[k].NumSparse());

		// as a sparse filter
		// convolution
		double	value_filter;
		for (int si=0; si<p_filters_flipped[k].NumSparse(); si++)
		{
			// site at flipped filter
			p_filters_flipped[k].GetValue(si,value_filter,j3,i3);

			j1 = half_row-j3; i1 = half_col-i3;

			// site at filtering result
			j2 = j+j1; i2 = i+i1;

			vc_resp(si) = p_filtering_result[k](j2,i2)-value_filter*(ras_image(j,i)-value);
		}
	}
	else
	{	// as a general filter
		// how many response changes regarding to the change at this site
		vc_resp.SetRowVector((row2-row1+1)*(col2-col1+1));

		// convolution
		int t = 0;
		for (j1=row1; j1<=row2; j1++)
			for (i1=col1; i1<=col2; i1++)
			{
				// site at filtering result
				j2 = j+j1; i2 = i+i1;
				// site at flipped filter
				j3 = half_row-j1; i3 = half_col-i1;

				vc_resp(t) = p_filtering_result[k](j2,i2)-p_filters_flipped[k](j3,i3)*(ras_image(j,i)-value);
				t++;
			}
	}
}

inline void MCMC_FilterBank::NewFilterRespAtOneSite(Histogram &histogram, const int k, const int j, const int i, const double value)
{
	int		j1,i1,j2,i2,j3,i3;
	int		half_row,half_col,row1,row2,col1,col2;
	double	temp;

	half_row = (p_filters_flipped[k].rows()-1)/2;
	half_col = (p_filters_flipped[k].cols()-1)/2;

	row1 = max(-half_row,0-j); row2 = min(half_row,ras_image.rows()-1-j);
	col1 = max(-half_col,0-i); col2 = min(half_col,ras_image.cols()-1-i);

	if (p_filters_flipped[k].NumSparse()<NUM_NONZERO_BINS)
	{
		double	value_filter;
		for (int si=0; si<p_filters_flipped[k].NumSparse(); si++)
		{
			// site at flipped filter
			p_filters_flipped[k].GetValue(si,value_filter,j3,i3);

			j1 = half_row-j3; i1 = half_col-i3;

			// site at filtering result
			j2 = j+j1; i2 = i+i1;

			temp = p_filtering_result[k](j2,i2)-value_filter*(ras_image(j,i)-value);
			histogram.AddBin(temp,1.0);
		}
	}
	else
	{
		// convolution
		for (j1=row1; j1<=row2; j1++)
			for (i1=col1; i1<=col2; i1++)
			{
				// site at filtering result
				j2 = j+j1; i2 = i+i1;
				// site at flipped filter
				j3 = half_row-j1; i3 = half_col-i1;

				temp = p_filtering_result[k](j2,i2)-p_filters_flipped[k](j3,i3)*(ras_image(j,i)-value);
				histogram.AddBin(temp,1.0);
			}
	}
}

inline void MCMC_FilterBank::UpdateFilterRespAtOneSite(const int j, const int i, const double value)
{
	int				j1,i1,j2,i2,j3,i3,k;
	int				half_row,half_col,row1,row2,col1,col2;

	for (k=0; k<num_filter; k++)
	{
		half_row = (p_filters_flipped[k].rows()-1)/2;
		half_col = (p_filters_flipped[k].cols()-1)/2;

		row1 = max(-half_row,0-j); row2 = min(half_row,ras_image.rows()-1-j);
		col1 = max(-half_col,0-i); col2 = min(half_col,ras_image.cols()-1-i);

		if (p_filters_flipped[k].NumSparse()<NUM_NONZERO_BINS)
		{
			double	value_filter;
			for (int si=0; si<p_filters_flipped[k].NumSparse(); si++)
			{
				// site at flipped filter
				p_filters_flipped[k].GetValue(si,value_filter,j3,i3);
				j1 = half_row-j3; i1 = half_col-i3;

				// site at filtering result
				j2 = j+j1; i2 = i+i1;

				p_filtering_result[k](j2,i2) = p_filtering_result[k](j2,i2)-value_filter*(ras_image(j,i)-value);
			}
		}
		else
		{
			// convolution
			for (j1=row1; j1<=row2; j1++)
				for (i1=col1; i1<=col2; i1++)
				{
					// site at filtering result
					j2 = j+j1; i2 = i+i1;
					// site at flipped filter
					j3 = half_row-j1; i3 = half_col-i1;

					p_filtering_result[k](j2,i2) = p_filtering_result[k](j2,i2)-p_filters_flipped[k](j3,i3)*(ras_image(j,i)-value);
				}
		}

	}
	ras_image(j,i) = value;
}


#endif
